home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlaev2.f < prev    next >
Text File  |  1997-06-25  |  5KB  |  171 lines
  1.       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
  10. *     ..
  11. *
  12. *  Purpose
  13. *  =======
  14. *
  15. *  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
  16. *     [  A   B  ]
  17. *     [  B   C  ].
  18. *  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
  19. *  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
  20. *  eigenvector for RT1, giving the decomposition
  21. *
  22. *     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
  23. *     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
  24. *
  25. *  Arguments
  26. *  =========
  27. *
  28. *  A       (input) DOUBLE PRECISION
  29. *          The (1,1) element of the 2-by-2 matrix.
  30. *
  31. *  B       (input) DOUBLE PRECISION
  32. *          The (1,2) element and the conjugate of the (2,1) element of
  33. *          the 2-by-2 matrix.
  34. *
  35. *  C       (input) DOUBLE PRECISION
  36. *          The (2,2) element of the 2-by-2 matrix.
  37. *
  38. *  RT1     (output) DOUBLE PRECISION
  39. *          The eigenvalue of larger absolute value.
  40. *
  41. *  RT2     (output) DOUBLE PRECISION
  42. *          The eigenvalue of smaller absolute value.
  43. *
  44. *  CS1     (output) DOUBLE PRECISION
  45. *  SN1     (output) DOUBLE PRECISION
  46. *          The vector (CS1, SN1) is a unit right eigenvector for RT1.
  47. *
  48. *  Further Details
  49. *  ===============
  50. *
  51. *  RT1 is accurate to a few ulps barring over/underflow.
  52. *
  53. *  RT2 may be inaccurate if there is massive cancellation in the
  54. *  determinant A*C-B*B; higher precision or correctly rounded or
  55. *  correctly truncated arithmetic would be needed to compute RT2
  56. *  accurately in all cases.
  57. *
  58. *  CS1 and SN1 are accurate to a few ulps barring over/underflow.
  59. *
  60. *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  61. *  Underflow is harmless if the input data is 0 or exceeds
  62. *     underflow_threshold / macheps.
  63. *
  64. * =====================================================================
  65. *
  66. *     .. Parameters ..
  67.       DOUBLE PRECISION   ONE
  68.       PARAMETER          ( ONE = 1.0D0 )
  69.       DOUBLE PRECISION   TWO
  70.       PARAMETER          ( TWO = 2.0D0 )
  71.       DOUBLE PRECISION   ZERO
  72.       PARAMETER          ( ZERO = 0.0D0 )
  73.       DOUBLE PRECISION   HALF
  74.       PARAMETER          ( HALF = 0.5D0 )
  75. *     ..
  76. *     .. Local Scalars ..
  77.       INTEGER            SGN1, SGN2
  78.       DOUBLE PRECISION   AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
  79.      $                   TB, TN
  80. *     ..
  81. *     .. Intrinsic Functions ..
  82.       INTRINSIC          ABS, SQRT
  83. *     ..
  84. *     .. Executable Statements ..
  85. *
  86. *     Compute the eigenvalues
  87. *
  88.       SM = A + C
  89.       DF = A - C
  90.       ADF = ABS( DF )
  91.       TB = B + B
  92.       AB = ABS( TB )
  93.       IF( ABS( A ).GT.ABS( C ) ) THEN
  94.          ACMX = A
  95.          ACMN = C
  96.       ELSE
  97.          ACMX = C
  98.          ACMN = A
  99.       END IF
  100.       IF( ADF.GT.AB ) THEN
  101.          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
  102.       ELSE IF( ADF.LT.AB ) THEN
  103.          RT = AB*SQRT( ONE+( ADF / AB )**2 )
  104.       ELSE
  105. *
  106. *        Includes case AB=ADF=0
  107. *
  108.          RT = AB*SQRT( TWO )
  109.       END IF
  110.       IF( SM.LT.ZERO ) THEN
  111.          RT1 = HALF*( SM-RT )
  112.          SGN1 = -1
  113. *
  114. *        Order of execution important.
  115. *        To get fully accurate smaller eigenvalue,
  116. *        next line needs to be executed in higher precision.
  117. *
  118.          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  119.       ELSE IF( SM.GT.ZERO ) THEN
  120.          RT1 = HALF*( SM+RT )
  121.          SGN1 = 1
  122. *
  123. *        Order of execution important.
  124. *        To get fully accurate smaller eigenvalue,
  125. *        next line needs to be executed in higher precision.
  126. *
  127.          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  128.       ELSE
  129. *
  130. *        Includes case RT1 = RT2 = 0
  131. *
  132.          RT1 = HALF*RT
  133.          RT2 = -HALF*RT
  134.          SGN1 = 1
  135.       END IF
  136. *
  137. *     Compute the eigenvector
  138. *
  139.       IF( DF.GE.ZERO ) THEN
  140.          CS = DF + RT
  141.          SGN2 = 1
  142.       ELSE
  143.          CS = DF - RT
  144.          SGN2 = -1
  145.       END IF
  146.       ACS = ABS( CS )
  147.       IF( ACS.GT.AB ) THEN
  148.          CT = -TB / CS
  149.          SN1 = ONE / SQRT( ONE+CT*CT )
  150.          CS1 = CT*SN1
  151.       ELSE
  152.          IF( AB.EQ.ZERO ) THEN
  153.             CS1 = ONE
  154.             SN1 = ZERO
  155.          ELSE
  156.             TN = -CS / TB
  157.             CS1 = ONE / SQRT( ONE+TN*TN )
  158.             SN1 = TN*CS1
  159.          END IF
  160.       END IF
  161.       IF( SGN1.EQ.SGN2 ) THEN
  162.          TN = CS1
  163.          CS1 = -SN1
  164.          SN1 = TN
  165.       END IF
  166.       RETURN
  167. *
  168. *     End of DLAEV2
  169. *
  170.       END
  171.